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SIGNAL PROCESSING 

Tlie present invention relates to the reconstruction of 
signals. There are many applications, such as radar, 
5 sonar, acoustic data, spectroscopy, geophysical, and image 
signal processing, in which it is desirable to reconstruct 
signals from given data. 

In many of these situations the effect on the signals 
of noise and the particular characteristic of the system 
10 generating the signal are known or can be approximated 
using appropriate mathematical models* In these 

•situations, Bayesian reconstruction methods have often been 
applied to reconstruct the signal from given data. These 
methods can work well. For example one Bayesian 

15 reconstruction approach, known as the maximum entropy 
method, is known to work well. Usually, the maximum 
entropy method (MEM) , is only be applied to the 
reconstruction of signals that are strictly non-negative, 
and which are not complex. Given this, it is generally not 
20 possible to change the basis in any vector space 
representing the data during any reconstruction. Being 
unable to do this results in computationally very expensive 
reconstruction processes having to be employed, because of 
the requirement for very large calculations having an 
25 extremely large nxamber of variables to be determined. 

For example, employment of such an MEM to a stack of 
twenty images from a microscope takes in the region of 

TM TH 

fifty minutes using a standard INTEL Pentium 200MHz 
processor. 

3 0 The present invention is directed towards improving 

the reconstruction of signals from given data so that such 
reconstruction can be performed within a reduced time 
frame . 

According to the present invention there is provided 
35 a method of reconstructing a signal from a given set of 
data, with a prediction function representing a predictable 



2 



effect on the signal, and a noise function representing 
unpredictable noise, the method comprising the steps of: 

altering the coordinate basis of the data and signal 
from an original coordinate basis in order to produce a 
prediction function having a reduced set of variables; 

performing a Bayesian reconstruction capable of 
operation of positive, negative, and complex signal values 
to produce a reconstruction signal; and 

converting the reconstruction signal back into the 
original coordinate basis to generate a signal. 

The Bayesian reconstruction may be performed using a 
Fourier basis, or may use a wavelet basis. 

The Bayesian reconstruction may employ the maximiam 
entropy method, in which case the method may employ an 
evaluation parameter, oc, which may be deteirmined from a 
prior reconstruction, set at a fixed value, or determined 
during the reconstruction step- 

The signal to be reconstructed may be an image signal, 
or may be a radar, sonar, or acoustic data signal. 
Alternatively, it may be a signal from spectroscopy or a 
geophysical signal. 

By employing the method of the present invention, an 
example stack of 2 0 microscope images takes approximately 
45 seconds to reconstruct using a standard INTEL Pentium 
200 MH^ processor. 

One example of the present invention will now be 
described. 
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Bayesian reconstruction methods have been applied to numerous problems in a wide variety of fields. In 
their standard form, however, they can be very computationally intensive, since they generally require the 
numerical maximisation of a complicated function of many variables. For example, in image reconstruction 
problems, it is not unusual for the number of variables to be ^ 10^. Furthermore, one of the most popular 
Bayesian reconstruction algorithms is the maximum-entropy method (MEM), which can only be applied 
to the reconstruction of signals that are strictly non-negative (see below). The method can, however, be 
extended to signals that can take both positive and negative values. We develop the MEM approach -so that 
it can be applied to the reconstruction of signals that can take positive, negative or complex values. As a 
result, this enables the use of similarity transformations in the reconstruction algorithm so that calculations 
can be performed in an alternative *basis' that is more appropriate to the problem under consideration. 
Specifically, the basis is chosen so that signal is reconstructed by performing a large number of numerical 
maximisations of low dimensionality, rather than a single maximisation of high dimensionality. This results 
in a significant increase in speed. Indeed, in the example outlined below, the speed of the reconstruction 
algorithm is increased by a factor of about 100. 

»: standard Bayesian reconstruction techniques mentioned above are are typically applied to a given data 
d{y) in order to reconstruct some underlying signal ^(x). Here, y denotes the space over which the data 
are defined and x denotes the space of the signal, which may in general be distinct from x and need not have 
the same number of dimensions. For example, given some data stream d(t), which consists of the measured 
values of some quantity as a function of time <, we may wish to reconstruct the two-dimensional spatial 
variation of some other quantity (or signal) s{x, y). 

For most measurements it is convenient (or necessary) to digitise the data and the signal (for example if 
either is to be stored/analysed on a digital computer). We may therefore denote the data by the vector d 
with Nti components, where Nd is the number of data samples. Similarly, we denote the signal by the vector 
s of length Ng , where Ns is the number of points at which we wish to reconstruct the signal. 

In general, we may express the data vector as some function © of the signal vector, i.e. 

d = e(s). 

The function 0 can be non-linear and specifies the effect of the measuring apparatus on the signal that 
we wish to reconstruct. It is customary to divide this function into the predictable effect of the measuring 
apparatus on the signal and the stochastic noise part due to inherent inaccuracies in the measurement 
process. In this case, we may write 

d = ^(s) + e, (1) 



wt»v=A e ^ denotes the predictable res^BR of the apparatus to the signal and e is ^^Ktor of length Nd that 
contains any stochastic noise contributions to the data. 

The Bayesian approach to reconstructing the signal is to calculate the estimator s that maximises the 
posterior probability Pr(s|d). This is given by Bayes^ theorem as 

Pr(s|d)- p^^jj , 

where Pr(d|s) is the likelihood of obtaining the data given the signal, Pr(s) is the prior probability, and the 
evidence Pr(d) can be considered merely as a normalisation constant. Thus, in order to obtain the Bayesian 
estimator of the signal vector, we must maximise the product Pr(d|s) Pr(s) of the likelihood function and 
the prior. A thorough discussion of Bayesian analysis techniques is given by Sivia (1996). 

The likelihood function describes the statistics of the noise contribution e to the data. This function may 
take any form appropriate to the noise statistics. It is convenient to define the log-likelihood function 
L{s) = ln[Pr(d|s)] so that the likelihood function may be written as Pr(d|s) = exp[Z/(s)]. As an example, 
if the noise on the data is Gaussian-distributed and described by the noise covariance matrix N , then the 
likelihood function takes the form 

Pr(d|s) oc exp[-|cTN-ic] 

oc exp [-l(d - ^(s))'rN-i (<* - n^))] . (2) 

where in the second line we have used (1). In this case, the log-likelihood function is simply minus one-half 
of the standard misfit statistic, i.e. L{s) = — |x^(s)- 

The prior distribution Pr(s) codifies our knowledge of the underlying signal before acquiring the data. If we 
have some advance knowledge of the statistical properties of the signal then this is contained in the prior. 
For example, if we assume the signal to be described by a Gaussian random field with a covariance matrix 
C, then the prior takes the form 

Pr(s) cc exp(-|sTC-^s). 

Indeed, if the prior is assumed to have this form and the likelihood is also Gaussian, as in (2), then the 
Bayesian estimator s obtained by maximising their product is identical to the standard Wieijer filter solution. 
An introduction to the Wiener filter technique is given by Press et al. (1994) 

It is clear, however, that although the noise contribution c to the data may often be Gaussian-distributed, 
the assumption of a Gaussian form for the prior is not valid for a general signal. If the joint probability 
distribution of the elements of the signal vector is known then it should be used as the prior. This is almost 
always impossible, however, and we instead investigate the assignment of a prior applicable to general 
signals that is based on information-theoretic considerations alone. Using very general notions of subset 
independence, coordinate in variance and system independence, it may be shown that the prior probability 
Pr(tf) should take the form 

Pr(s) oc exp[a5(s, m)], (3) 

where the dimensional constant a depends on the scaling of the problem and may be considered as a 
regularising parameter, and m is a model to which the Bayesian reconstruction defaults in the absence of 
any data and is usually set to a small constant value. The function 5(s, m) is the cross-entropy of the signal 
and model vectors and is given by 

5(s, m) = ^ [*„ - mn - fi„ In (^)] , (4) 

where AT, is the length of the (digitised) signal vector. A derivation of this result is given by Skilling (1988). 
By combining the entropic prior with the likelihood function, the Bayesian estimator of the signal is found 
by maximising with respect to s the posterior probability, which now takes the form 

Pr(s|d) oc exp[i.(s)] exp[a5(s, m)] = exp[L(s) 4- cr5(s, m)]. 

Thus, maximising this probability distribution is equivalent to maximising the function 

P(s) = L(s) + a5(s,m), (5) 



auv* chis forms the basis of the maximum entropy method (MEM). 

The maximum-entropy method has been applied to a wide range of signal reconstruction problems. In 
its standard form, however, it can be very computationally intensive. The function F{s) is in general a 
complicated function of the components Sn of the signal vector and so a numerical maximisation of the F{s) 
must be perform over this iV, -dimensional space. It is not unusual for to be of the order ^ 10®, 
particularly in image reconstruction problems. Moreover, the standard MEM approach is only applicable to 
signals that are strictly non-negative, as is clear from the presence of the logarithmic term in the expression 
(4) for the entropy. 

Nevertheless, it is possible to extend the MEM so that it can be applied to the reconstruction of signals 
that can take both positive and negative values. The definition of the entropy for positive/negative signals 
with certain special properties was first presented by Gull & Skilling (1990). The generalistion to arbitrary 
positive/negative signals and the derivation of the prior probability in this case is given by Hobson & Lasenby 
(1998). It is found that the prior has the same form as given (3) but the expression for the entropy 5(s, m) 
must be modified. The central idea is to express the general (positive/negative) signal vector s as the 
difference of two vectors u and v containing only strictly non-negative distributions, i.e. 



s = u — v. 
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applying continuity constraints on the entropy functional it is possible to show that the expression for 
entropy for the positive/negative signal s is given by 

S(s, m„, m„) = 53 - {rnu)n - (mv)„ - s„ In | » (^) 

where mu and m„ are separate models for u and v respectively, and where Vn = + 4(mu)n(mv)„]^/^. We 
cannot hope to replace the models m„ and mt, by a single positive/negative model m, (say), since such a 
distribution could no longer be considered as an integration measure. Nevertheless, we can still consider the 
difference m« - m,, as the model for the signal s. We note that the form of the positive/negative entropy 
derived by Gull & Skilling (1990) requires m„ = mt,. 

Given the entropic prior for general positive/negative signals, it is then straightforward to define the prior 
for complex signals simply by applying the above analysis to the real and imaginary parts separately. In 
this case the models m,^ and mv are also taken to be complex. The real and imaginary parts of are the 
models for the positive portions of the real and imaginary parts of s respectively. Similarly, the real and 
imaginary parts of are the models for the negative portions of the real and imaginsiry parts of the image. 
The total entropy of the complex signal is then obtained by evaluating the sum (6) using first the real parts 
and then the imaginary parts of s, nriu and mv, and adding the results. 

The ability to reconstruct positive/negative and complex distributions using the MEM approach has pro- 

«md consequences for greatly improving the both the speed and accuracy of the MEM technique. These 
provements are based on the idea of making a change of basis in both the signal and data spaces and 
performing a Bayesian reconstruction in this new basis. With an appropriate choice of the new bases, it 
is possible to speed up significantly the calculation of the reconstruction, which can then easily be rotated 
back into the original basis in signal space to obtain the final reconstruction of the signal. For example, 
the ability to reconstruct complex signals allows us to perform reconstructions in the Fourier basis of com- 
plex exponentials, which greatly reduces the computational complexity of de-blurring images that have been 
convolved with a spatially varying point-spread function (see the example below). 

In order to understand how a general change of basis is performed we must first remind ourselves of some 
basic results in linear algebra and vector spaces (see e.g. Riley, Hobson & Bence 1997). Suppose there exists 
a set of linearly independent vectors e^") (n = 1, . . iV^,) that form a complete basis for the iSA, -dimensional 
space of the signal vector. We may then write the signal vector as a weighted sum of these basis vectors. 
Formally, if we take e^") to be the column vector with unity as the nth element and zeros elsewhere, then 
we may write the signal vector as 



= E 



(n) 



.fine- 

Thus we see that in order to reconstruct the signal vector, we are in fact reconstructing its coefficients in 
this trivial basis. We can, however, equally well expand the signal vector in terms of any other linearly 



inuependent set of vectors (n = I^^, N,) such that 



We may perform a similar procedure in the JSr^-dimensional data space, which is in general distinct from the 
signal space. If we consider the trivial basis vectors f (n = 1, . . . , Nd) in this space, with unity as the nth 
elements and zeros elsewhere, then the data vector is given by 

On performing change of basis in data space to some other basis f this becomes 

Since the noise vector e also belongs to the data space a similar change of basis can applied to it, such that 

n=l n=l 

It is clear that, even if the elements Sn of the signal vector in the original basis were strictly non-negative, 
the elements s'^ in the new basis will in general take both positive and negative values. Furthermore, in the 
case where the new basis vectors have complex components, the coefficients s'^ may themselves be complex. 
Hence it is the extension of the MEM tedinique to the reconstruction of such quantities that allows this 
approach to be taken. 

Once we have performed the changes of basis in the data and signal spaces, we denote the vectors with 
components s'^ by s' and we similarly define the vectors d' and e' as those containing the elements and 
cjj respectively. In the signal space we can relate the two bases e(") and e'<") (n = 1, . . . , AT,) by 

where U{n is the ith component of e'^"^ with respect to the original (unprimed) basis. The vectors s' and s 
axe then related by 

s = Us'. (7) 
Similar results hold for the bases f and f'<"^ (n = 1, . . , , Nd) in data space, such that 

d = Vd'. (8) 

where the element Kn is the tth component of f'^"^ with respect to the unprimed basis. A similar expression 
exists relating the noise vectors c' and €. Substituting (7) and (8) into (1), we then obtain 

d' = -*'(s')+e', (9) 
where ^' is a new function relating the signal and data vectors in the new basis. 

Clearly, our aim is to choose the new bases in the data and signal spaces in order that the relationship (9) takes 
the simplest possible form. More formally, we wish to perform similarity transformations in the data and 
signal spaces that partition each space into numerous (quasi-)disjoint subspaces of much lower dimensionatlity. 
In such bases, we may then calculate the estimate s' of the rotated signal vector by performing numerical 
maximisations in each subspace separately. Thus, we replace a single maximisation over the iVj -dimensional 
signal space in the original basis by numerous low-dimensionalty maximisations over each subspace in the 
new basis. This leads to a considerable increase in the speed with which the reconstruction can be performed. 
Then, having calculated the Bayesian reconstruction s' in the new basis, the required signal reconstruction 
s can be obtained by rotating back to the original signal-space basis. We reiterate that, since the elements 
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s„ m the new basis can in general take positive, negative or even complex values, it is the extension of the 
MEM technique to the reconstruction of such quantities that allows this approach to be taken. 

As an example, let us consider the application of the above Bayesian reconstruction technique to the decon- 
volution of images that have been blurred by a spatially-invariauit point-spread function (PSF) and that may 
also contain some noise contribution. For simplicity, we will assume that the de-blurred reconstruction is 
produced on the same pixelisation as the blurred image, although this is clearly not required by the technique 
in general. In this simple case, the data and signal spaces coincide. 

It is well-known that the convolution of an underlying image with a spatially-invariant PSF is equivalent to 
multiplying together their Fourier transforms and then performing an inverse Fourier transform. Therefore, 
in the Fourier domain, each Fourier mode can be considered independently of the others. This suggests that 
we should perform the Bayesian reconstruction in the Fourier basis, such that 

«ft = 5Z exp(-t7rn(m - l)/7Nr,)sm. 

with a similar expression relating the components of the data vectors d' and d and the noise vectors c and e' 
(since the data and signal spaces coincide). Thus, in this case, the JV^ -dimensional signal (and data) space 
has been partitioned into Ng separate disjoint spaces (i.e. one for each Fourier mode), 

tfor each value of n (or Fourier mode), we may consider the elements dj,, s'^ and independently of 
; for other values of n. This leads to a substantial decrease in the CPU time required to de-blur a given 
image. For simplicity, at our chosen Fourier mode we denote d'^ s'„ and cj^ by d', s' and c' respectively. The 
quantity is given simply by the Fourier coefficient of the true underlying image, or signal s', mutliplied by 
the Fourier coefficient of the PSF, or response R, In addition, a noise contribution, c*, in the Fourier domain 
may also be present. If is no instrumental errors are expected from a given apparatus, it is still possible to 
introduce *noise* by, for example, digitising an image in order to store it on a computer. Thus the data value 
is given by 

rf' = -h c', (10) 

Since we eire performing the reconstruction in the Fourier basis, the noise on each Fourier mode will contain 
contributions from a wide range of scales. Therefore, provided the noise on the image is distributed in a 
statistically-homogeneous manner, we would expect from the central limit theorem that the noise in the 
Fourier domain is described reasonably well by a Gaussian distribution. Therefore, the likelihood function 
is given by 

Pr(d'|sO oc exp (-[d' - Rs'\^/a'') , 
where = W^*) ^ the variaince of the noise contribution at the particular Fourier mode under consideration. 
From (6), the entropy 5(^,m) of this complex 'image' may be shown to be given by (where we have set 

S(s, m) = - 25R(m) - In + - 2a(m) - S>(.') In [ ^^^^jm/^^ ] ' 

where the SR and O denote the real or imaginary part respectively of a complex number; also ^{if) = 
[3f?(s')^ -l-45R(m)2]^/^ ^ similar expression exists for 5>(^). 

Using the above expressions for the likelihood and prior, and assuming a particular value for the regularising 
parameter a in (5), it is then possible numerically to maximise the posterior probability to obtain the 
estimator P of the signal vector at each Fourier mode independently. Once these estimators have been 
calculated for all the Fourier modes, we simply perform an inverse Fourier transform to recover the de- 
blurred image, i.e. 



= }r S exp(+i7rn(m - 1)/N,)s'^ 



The value of a used in the reconstruction algorithm may be set in three different ways. Firstly, a may be set 
such that the misfit statistic between the observed data and that predicted from the reconstruction, is 
equal to its expectation value, i.e. the number of data points to be fitted. This choice is usually referred to 
as 'historic' MEM. Alternatively, it is possible to determine the appropriate value for a in a fully Bayesian 
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manner (Skilling 1989) by simply treating it as another parameter in our hypothesis space. This is the 
recommended choice and is usually referred to as ^classic* MEM. Finally, the simplest option is to iix a in 
advance to some value. This option is unlikely to yield optimal results. It can, however, be used to obtain a 
quick solution if the historic or classic value for a has already been determined for a particular problem on 
a previous occasion. 

The above technique has been applied to several different data 

sets in which an image has been convolved with a spatially-invariant point-spread function. Fajc, e^csnpXe 
the ^classic' . de-blurring of a microscope image of a section through a pollen grain. The section 

has dimensions 128 x 128 and is taken from a three-dimensional 'stack' of 20 such images. The original 
stack had been blurred by a spatially-invariant three-dimensional PSF and the de-blurred reconstruction of 
the entire stack required approximately 45 seconds on an Intel Pentium 200 MHz processor. A standard 
MEM algorithm, which does not use the similarity transformation technique, was also applied to this stack of 
images. This produced reconstructions of a similar quality to those obtained using tte inwenticn, but required 
approximately 50 minutes CPU time on the same machine. 

Similar gains in speed can be obtained by expanding the signal in different bases appropriate to the given 
problem, so long as the correlation matrix of the resulting coefficients of the data, signal and noise vectors 
in the new bases is relatively sparse. For example, similar results may be obtained by reconstructing the 
coefficients in a wavelet expansion (Daubechie 1992) of the signal as opposed to the Fourier expansion used 
above. This case has the advantage that the coefficients are always real. Furthermore, the scaling/translation 
properties of the wavelet transform allow automatic multi-resolution reconstructions of the signal vector. 

Clearly, the general method outlined above can be applied to numerous different reconstruction problems of 
arbitrary dimensionality. Examples include the analysis of acoustic data, radar, underwater sonar, spectro- 
scopy, geophysical data, oil exploration and medical imaging. In addition to spatial dimensions, additional 
dimensions such as time, spectral behaviour and polarisation are also easily accommodated. 
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CLAIMS 

1. A method of reconstructing a signal from a given set 
of data, with a prediction function representing a 

5 predictable effect on the signal, and a noise function 
representing unpredictable noise, the method comprising the 
steps of : 

altering the coordinate basis of the data and signal 
from an original coordinate basis in order to produce a 
10 prediction function having a reduced set of variables; 

performing a Bayesian reconstruction capable of 
operation of positive, negative, and complex signal values 
to produce a reconstruction signal; and 

converting the reconstruction signal back into the 
15 original coordinate basis to generate a signal. 

2. A method according to claim 1, wherein the Bayesian 
reconstruction is performed using a Fourier basis. 

20 3. A method according to claim 1, wherein the Bayesian 
reconstruction is performed using a wavelet basis. 

4. A method according to any preceding claim, wherein the 
Bayesian reconstruction employs the maximum entropy method. 



25 
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5. A method according to claim 4, employing an evaluation 
parameter, oc, which is determined from a prior 
reconstruction . 

6. A method according to claim 4, employing an evaluation 
parameter, oc, which is set at a fixed value. 



7. A method according to claim 4, employing an evaluation 
parameter, oc, which is deteirmined during the reconstruction 
35 step. 
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8. A method according to any of claims 1 to 7, in which 
the signal to be reconstructed is an image signal. 

9. A method according to claim 8, wherein the image 
5 signal is a medical image signal. 

10. A method according to any of claims 1 to 7, in which 
the signal to be reconstructed is a radar signal. 

10 11. A method according to any of claims 1 to 7, in which 
the signal to be reconstructed is an acoustic data signal . 

12. A method according to claim 11, wherein the acoustic 
data signal is an underwater sonar signal. 

15 

13. A method according to claim 11, wherein the acoustic 
data signal is a geophysical data signal. 

14. A method according to any of claims 1 to 7, in which 
20 the signal to be reconstructed is a signal from 

spectroscopy . 

15. A method according to any one of claims 1 to 7, in 
which the signal is a communication signal, such as a 

2 5 time-series signal. 
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Abstract 

A method of reconstructing a signal from a given set 
of data, with a prediction function representing a 
5 predictable effect on the signal, and a noise function 
representing unpredictable noise. The method comprises 
the steps of altering the coordinate basis of the data and 
signal from an original coordinate basis in order to 
produce a prediction function having a reduced set of 
10 variables, performing a Bayesian reconstruction capable of 
operation of positive, negative, and complex signal values 
to produce a reconstruction signal, and converting the 
reconstruction signal back into the original coordinate 
basis to generate a signal. 

15 



t 




t 




